Model annual dissimilarity with temperature, fishing, and survey identity

Script 8 for Kitchel et al. 2024 in prep taxonomic diversity manuscript.

library(data.table)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
data.table 1.14.8 using 1 threads (see ?getDTthreads).  Latest news: r-datatable.com
**********
This installation of data.table has not detected OpenMP support. It should still work but in single-threaded mode.
This is a Mac. Please read https://mac.r-project.org/openmp/. Please engage with Apple and ask them for support. Check r-datatable.com for updates, and our Mac instructions here: https://github.com/Rdatatable/data.table/wiki/Installation. After several years of many reports of installation problems on Mac, it's time to gingerly point out that there have been no similar problems on Windows or Linux.
**********

library(MuMIn)
library(ggplot2)
library(cowplot)
library(lme4)

Loading required package: Matrix

library(stringr)
library(nlme)

Attaching package: ‘nlme’

The following object is masked from ‘package:lme4’:

    lmList

#pull in function to calculate model estimates and standard errors
source(here::here("analysis_code","extract_coefficients_function.R"))

###Predicts annual dissimilarity with annual characteristics, temperature and fishing values

Pull in - region areas (if not already loaded) - region characteristics (if not already loaded); saveRDS(FishGlob_richness_year_survey, file = here::here(“output”,“FishGlob_richness_year_survey.Rds”)) - fishing (if not already loaded) - temp (if not already loaded)

#physical area by year
region_area_byyear <- fread(here::here("output","region_area_byyear.csv"))

#merged fishing, temp, dissimilarities
dissimilarities_temp_fishing <- fread(here::here("output","dissimilarities_temp_fishing.csv"))

#combine
dissimilarities_temp_fishing_area <- dissimilarities_temp_fishing[region_area_byyear, on = c("survey_unit","year")]

#only jaccard for these analyses
dissimilarities_temp_fishing_area.jaccard <- dissimilarities_temp_fishing_area[dissimilarity_metric == "jaccard_dissimilarity_index_binary",]

Add in season

#load up julian days
dates_regions <- readRDS(here::here("output","dates_regions.rds"))

#most common season per year
# Function to get the most frequent value
get_mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

# Apply the function to each group
dates_regions.r <- dates_regions[, .(season = get_mode(season)), by = c("survey_unit","year")]


#merge
dissimilarities_temp_fishing_area.jaccard <- dates_regions.r[dissimilarities_temp_fishing_area.jaccard, on = c("survey_unit","year")]

Pull in palette and name helper

source(here::here("analysis_code","color_links.R"))

Pull in observed trend values

jaccard_total_coefs.r <- fread(here::here("output","jaccard_total_coefs.r.csv"))

Plot fishing and temperature vs. time for all regions

#######TEMPERATURE

#set order by survey unit for plotting
all_surveys <- levels(as.factor(dissimilarities_temp_fishing_area.jaccard$survey_unit))
setorder(dissimilarities_temp_fishing_area.jaccard, survey_unit)

dissimilarities_temp_fishing_area.jaccard[,Survey_Name_Season:=factor(Survey_Name_Season, levels = unique(dissimilarities_temp_fishing_area.jaccard$Survey_Name_Season), ordered = T)]

(sbt_time_survey_facet_1_20 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard[survey_unit %in% all_surveys[1:20] & year > 1979]) +
  labs(y = "Mean bottom temperature (˚C)",  x = "Year") +
  geom_point(aes(y = as.numeric(yearly_mean_bypoint_avg), x = year), alpha = 0.3) +
  geom_smooth(aes(y = as.numeric(yearly_mean_bypoint_avg), x = year), method = "lm") +
  scale_x_continuous(breaks = ~ axisTicks(., log = FALSE)) +
  theme_classic() +
  theme(axis.text.x = element_text(size = 7)) +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4))

ggsave(sbt_time_survey_facet_1_20, path = here::here("figures"), filename = "sbt_time_survey_facet_1_20.jpg", height = 12, width =9)

(sbt_time_survey_facet_21_34 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard[survey_unit %in% all_surveys[21:34] & year > 1979]) +
  labs(y = "Mean bottom temperature (˚C)",  x = "Year") +
  geom_point(aes(y = as.numeric(yearly_mean_bypoint_avg), x = year), alpha = 0.3) +
  geom_smooth(aes(y = as.numeric(yearly_mean_bypoint_avg), x = year), method = "lm") +
  scale_x_continuous(breaks = ~ axisTicks(., log = FALSE)) +
  theme_classic() +
  theme(axis.text.x = element_text(size = 7)) +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4))

ggsave(sbt_time_survey_facet_21_34, path = here::here("figures"), filename = "sbt_time_survey_facet_21_34.jpg", height = 12, width =9)

#######FISHING

dissimilarities_temp_fishing_area.jaccard.cc <- dissimilarities_temp_fishing_area.jaccard[complete.cases(dissimilarities_temp_fishing_area.jaccard[,summed_tonnes_scaled_byreg]),]

(fishing_time_survey_facet_1_20 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard.cc[survey_unit %in% all_surveys[1:20] & year > 1979]) +
  labs(y = "Relative fishing catch",  x = "Year") +
  geom_point(aes(y = summed_tonnes_scaled_byreg, x = year), alpha = 0.3) +
  geom_smooth(aes(y = summed_tonnes_scaled_byreg, x = year), method = "lm") +
  scale_x_continuous(breaks = ~ axisTicks(., log = FALSE)) +
  theme_classic() +
  theme(axis.text.x = element_text(size = 7)) +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4))

ggsave(fishing_time_survey_facet_1_20, path = here::here("figures"), filename = "fishing_time_survey_facet_1_20.jpg", height = 12, width =9)

(fishing_time_survey_facet_21_34 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard.cc[survey_unit %in% all_surveys[c(21:34)] & year > 1979]) +
  labs(y = "Relative fishing catch",  x = "Year") +
  geom_point(aes(y = summed_tonnes_scaled_byreg, x = year), alpha = 0.3) +
  geom_smooth(aes(y = summed_tonnes_scaled_byreg, x = year), method = "lm") +
  scale_x_continuous(breaks = ~ axisTicks(., log = FALSE)) +
  theme_classic() +
  theme(axis.text.x = element_text(size = 7)) +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4))

ggsave(fishing_time_survey_facet_21_34, path = here::here("figures"), filename = "fishing_time_survey_facet_21_34.jpg", height = 12, width =9)

NA
NA

Plot fishing and temperature vs. dissimilarity for all regions

#####MEAN TEMP
(preds_sbt_mean_temp_survey_facet_1_20 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard[survey_unit %in% all_surveys[1:20] & year > 1979]) +
  labs(x = "Mean bottom temperature (˚C)",  y = "β-diversity") +
  geom_point(aes(x = as.numeric(yearly_mean_bypoint_avg), y = annual_dissimilarity_value), alpha = 0.3) +
  geom_smooth(aes(x = as.numeric(yearly_mean_bypoint_avg), y = annual_dissimilarity_value), method = "lm") +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4) +
  theme_classic())

ggsave(preds_sbt_mean_temp_survey_facet_1_20, path = here::here("figures"), filename = "preds_sbt_mean_temp_survey_facet_1_20.jpg", height = 12, width =9)

(preds_sbt_mean_temp_survey_facet_21_34 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard[survey_unit %in% all_surveys[21:34] & year > 1979]) +
  labs(x = "Mean bottom temperature (˚C)",  y = "β-diversity") +
  geom_point(aes(x = as.numeric(yearly_mean_bypoint_avg), y = annual_dissimilarity_value), alpha = 0.3) +
  geom_smooth(aes(x = as.numeric(yearly_mean_bypoint_avg), y = annual_dissimilarity_value), method = "lm") +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4) +
  theme_classic())

ggsave(preds_sbt_mean_temp_survey_facet_21_34, path = here::here("figures"), filename = "preds_sbt_mean_temp_survey_facet_21_34.jpg", height = 12, width =9)

#####MINIMUM TEMP
(preds_sbt_min_temp_survey_facet_1_20 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard[survey_unit %in% all_surveys[1:20] & year > 1979]) +
  labs(x = "Minimum bottom temperature (˚C)",  y = "β-diversity") +
  geom_point(aes(x = as.numeric(yearly_min_bypoint_avg), y = annual_dissimilarity_value), alpha = 0.3) +
  geom_smooth(aes(x = as.numeric(yearly_min_bypoint_avg), y = annual_dissimilarity_value), method = "lm") +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4) +
  theme_classic())

ggsave(preds_sbt_min_temp_survey_facet_1_20, path = here::here("figures"), filename = "preds_sbt_min_temp_survey_facet_1_20.jpg", height = 12, width =9)

(preds_sbt_min_temp_survey_facet_21_34 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard[survey_unit %in% all_surveys[21:34] & year > 1979]) +
  labs(x = "Minimum bottom temperature (˚C)",  y = "β-diversity") +
  geom_point(aes(x = as.numeric(yearly_min_bypoint_avg), y = annual_dissimilarity_value), alpha = 0.3) +
  geom_smooth(aes(x = as.numeric(yearly_min_bypoint_avg), y = annual_dissimilarity_value), method = "lm") +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4) +
  theme_classic())

ggsave(preds_sbt_min_temp_survey_facet_21_34, path = here::here("figures"), filename = "preds_sbt_min_temp_survey_facet_21_34.jpg", height = 12, width =9)

#######FISHING

dissimilarities_temp_fishing_area.jaccard.cc <- dissimilarities_temp_fishing_area.jaccard[complete.cases(dissimilarities_temp_fishing_area.jaccard[,summed_tonnes_scaled_byreg]),]

(preds_sbt_mean_fishing_survey_facet_1_20 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard.cc[survey_unit %in% all_surveys[1:20] & year > 1979]) +
  labs(x = "Relative fishing catch",  y = "β-diversity") +
  geom_point(aes(x = summed_tonnes_scaled_byreg, y = annual_dissimilarity_value), alpha = 0.3) +
  geom_smooth(aes(x = summed_tonnes_scaled_byreg, y = annual_dissimilarity_value), method = "lm") +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4) +
  theme_classic())

ggsave(preds_sbt_mean_fishing_survey_facet_1_20, path = here::here("figures"), filename = "preds_sbt_mean_fishing_survey_facet_1_20.jpg", height = 12, width =9)

(preds_sbt_mean_fishing_survey_facet_21_34 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard.cc[survey_unit %in% all_surveys[c(21:34)] & year > 1979]) +
  labs(x = "Relative fishing catch",  y = "β-diversity") +
  geom_point(aes(x = summed_tonnes_scaled_byreg, y = annual_dissimilarity_value), alpha = 0.3) +
  geom_smooth(aes(x = summed_tonnes_scaled_byreg, y = annual_dissimilarity_value), method = "lm") +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4) +
  theme_classic())

ggsave(preds_sbt_mean_fishing_survey_facet_21_34, path = here::here("figures"), filename = "preds_sbt_mean_fishing_survey_facet_21_34.jpg", height = 12, width =9)

NA
NA

###Plot number of tows per year per region

tows_year <- unique(dissimilarities_temp_fishing_area.jaccard[,.(survey_unit, haul_id_count_annual, area_km, year, Survey_Name_Season)])

min_haul_density <- min(tows_year$haul_id_count_annual/tows_year$area_km)
max_haul_density <- max(tows_year$haul_id_count_annual/tows_year$area_km)

(tow_density_survey_facet_1_20 <- ggplot(data = tows_year[survey_unit %in% all_surveys[1:20]]) +
  labs(x = "Year",  y = expression("Tow density (tows per km"^2*")")) +
  geom_point(aes(x = year, y = haul_id_count_annual/area_km)) +
  ylim(c(min_haul_density-0.0001, max_haul_density+0.0001)) +
  facet_wrap(~Survey_Name_Season, ncol = 4) +
  theme_classic())

ggsave(tow_density_survey_facet_1_20, path = here::here("figures"), filename = "tow_density_survey_facet_1_20.jpg", height = 12, width =9)

(tow_density_survey_facet_21_34 <- ggplot(data = tows_year[survey_unit %in% all_surveys[21:34]]) +
  labs(x = "Year",  y = expression("Tow density (tows per km"^2*")")) +
  geom_point(aes(x = year, y = haul_id_count_annual/area_km)) +
  ylim(c(min_haul_density-0.0001, max_haul_density+0.0001)) +
  facet_wrap(~Survey_Name_Season, ncol = 4) +
  theme_classic())

ggsave(tow_density_survey_facet_21_34, path = here::here("figures"), filename = "tow_density_survey_facet_21_34.jpg", height = 12, width =9)

#minimum and maximum dissimilarity value
min_dissimilarity <- min(dissimilarities_temp_fishing_area.jaccard$annual_dissimilarity_value)
max_dissimilarity <- max(dissimilarities_temp_fishing_area.jaccard$annual_dissimilarity_value)

#how does tow density vary with dissimilarity?
tow_density_dissimilarity_1_20 <- 
  ggplot(data = dissimilarities_temp_fishing_area.jaccard[dissimilarity_metric == "jaccard_dissimilarity_index_binary" & survey_unit %in% all_surveys[1:20]]) +
  geom_point(aes(x = (haul_id_count_annual/area_km), y = annual_dissimilarity_value, color = year), size = 2) +
  viridis::scale_color_viridis() +
  facet_wrap(~Survey_Name_Season, ncol = 4, scales = "free") +
  labs(x = expression("Tow density (tows per km"^2*")"), y = "β-diversity", color = "Year") +
   # xlim(c(min_haul_density-0.0001, max_haul_density+0.0001)) +
   # ylim(c(min_dissimilarity-0.01,max_dissimilarity+0.01)) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

ggsave(tow_density_dissimilarity_1_20, path = here::here("figures"), filename = "tow_density_dissimilarity_1_20.jpg", height = 12, width =9)

tow_density_dissimilarity_21_34 <- 
  ggplot(data = dissimilarities_temp_fishing_area.jaccard[dissimilarity_metric == "jaccard_dissimilarity_index_binary" & survey_unit %in% all_surveys[21:34]]) +
  geom_point(aes(x = (haul_id_count_annual/area_km), y = annual_dissimilarity_value, color = year), size = 2) +
  viridis::scale_color_viridis() +
  facet_wrap(~Survey_Name_Season, ncol = 4, scales = "free") +
  labs(x = expression("Tow density (tows per km"^2*")"), y = "β-diversity", color = "Year") +
   # xlim(c(min_haul_density-0.0001, max_haul_density+0.0001)) +
   # ylim(c(min_dissimilarity-0.01,max_dissimilarity+0.01)) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

ggsave(tow_density_dissimilarity_21_34, path = here::here("figures"), filename = "tow_density_dissimilarity_21_34.jpg", height = 12, width =9)

#Plot all at once?
tow_density_dissimilarity <- 
  ggplot(data = dissimilarities_temp_fishing_area.jaccard[dissimilarity_metric == "jaccard_dissimilarity_index_binary"]) +
  geom_point(aes(x = (haul_id_count_annual/area_km), y = annual_dissimilarity_value, color = Survey_Name_Season), size = 2) +
  scale_color_manual(values = color_link$hex) +
  labs(x = expression("Tow density (tows per km"^2*")"), y = "β-diversity", color = "Year") +
   # xlim(c(min_haul_density-0.0001, max_haul_density+0.0001)) +
   # ylim(c(min_dissimilarity-0.01,max_dissimilarity+0.01)) +
  theme_classic()

#and by numbers?
result <- dissimilarities_temp_fishing_area.jaccard[, .(max_density = max(haul_id_count_annual/area_km, na.rm = TRUE), 
                 min_density = min(haul_id_count_annual/area_km, na.rm = TRUE)), 
             by = Survey_Name_Season]

result[,range:= max_density-min_density]

ggplot(data = result) +
  geom_point(aes(x = Survey_Name_Season, y = range)) +
  theme_classic() +
  labs(x = "",y = "Range of tow density values in tow/km^2") +
  theme(axis.text.x = element_text(angle = 45, hjust=1))

Dissimilarity by season of sampling

dissimilarities_temp_fishing_area.jaccard[,season := factor(season, levels = c("Spring","Summer","Fall","Winter"))]

ggplot(dissimilarities_temp_fishing_area.jaccard) +
  geom_boxplot(aes(x = season, y = annual_dissimilarity_value)) +
  theme_classic()

###Set up dredge to identify best performing models

options(na.action = "na.fail")

dissimilarity_covariates_dredge.dt <- dissimilarities_temp_fishing_area.jaccard[,.
                  (year, survey_unit,
                    yearly_mean_bypoint_avg, yearly_max_bypoint_avg, yearly_min_bypoint_avg,yearly_seas_bypoint_avg,
                    yearly_mean_bypoint_SD, yearly_max_bypoint_SD, yearly_min_bypoint_SD,yearly_seas_bypoint_SD,
                    yearly_mean_bypoint_avg.s, yearly_max_bypoint_avg.s, yearly_min_bypoint_avg.s,yearly_seas_bypoint_avg.s,
                    annual_dissimilarity_value,
                    haul_id_count_annual,
                    spp_count_annual, depth_annual_avg,
                    depth_annual_range, latitude_annual_avg,
                    latitude_annual_range, area_km, season, summed_tonnes_scaled_byreg)]

#merge in with colors for plotting predictions by survey
dissimilarity_covariates_dredge.dt <- color_link[dissimilarity_covariates_dredge.dt, on = "survey_unit"]

#If NA for any covariate, delete row
#View(dissimilarity_covariates_dredge.dt)
#Deleted:
  #Before 1980 and after 2019
  #Gulf of Saint Laurence South (no depth data)
  #No clear SAU match for Rockall Plateau
dissimilarity_covariates_dredge.dt <- dissimilarity_covariates_dredge.dt[complete.cases(dissimilarity_covariates_dredge.dt)]


dissimilarity_covariates_dredge.dt[, yearly_mean_bypoint_avg.scaledacrossall := scale(yearly_mean_bypoint_avg)][, yearly_mean_bypoint_SD.scaledacrossall := scale(yearly_mean_bypoint_SD)][, yearly_min_bypoint_avg.scaledacrossall := scale(yearly_min_bypoint_avg)][, yearly_max_bypoint_avg.scaledacrossall := scale(yearly_max_bypoint_avg)][, yearly_seas_bypoint_avg.scaledacrossall := scale(yearly_seas_bypoint_avg)][,haul_id_count_annual.scaledacrossall := scale(haul_id_count_annual)][,spp_count_annual.scaledacrossall := scale(spp_count_annual)][,depth_annual_avg.scaledacrossall := scale(depth_annual_avg)][,depth_annual_range.scaledacrossall := scale(depth_annual_range)] [,latitude_annual_avg.scaledacrossall := scale(latitude_annual_avg)][,latitude_annual_range.scaledacrossall := scale(latitude_annual_range)][,area_km.scaledacrossall := scale(area_km)]

###Full model

temp + survey_unit + fishing + area + latitude range + latitude average + depth range + depth average + spp count + # of hauls + season + AR for year

For temperature, we will look at: -mean (scaled across all regions) -max (scaled across all regions) -min (scaled across all regions) -seas (scaled across all regions) -SD

Comparing temp variables AND the presence/absence random effect for year

global_mod_mean_temp_lmer_surveyyear <- lmer(annual_dissimilarity_value ~ 
                             survey_unit:yearly_mean_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit:summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                             season + yearly_mean_bypoint_avg.scaledacrossall + summed_tonnes_scaled_byreg +
                               (1|survey_unit) + (1|year),
                             data = dissimilarity_covariates_dredge.dt)

global_mod_mean_temp_lmer_survey <- lmer(annual_dissimilarity_value ~ 
                             survey_unit:yearly_mean_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit:summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                             season + yearly_mean_bypoint_avg.scaledacrossall + summed_tonnes_scaled_byreg +
                               (1|survey_unit),
                             data = dissimilarity_covariates_dredge.dt)

global_mod_max_temp_lmer_surveyyear <- lmer(annual_dissimilarity_value ~ 
                             survey_unit:yearly_max_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit:summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                             season + yearly_mean_bypoint_avg.scaledacrossall + summed_tonnes_scaled_byreg +
                               (1|survey_unit) + (1|year),
                             data = dissimilarity_covariates_dredge.dt)

global_mod_max_temp_lmer_survey <- lmer(annual_dissimilarity_value ~ 
                             survey_unit:yearly_max_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit:summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                             season + yearly_mean_bypoint_avg.scaledacrossall + summed_tonnes_scaled_byreg +
                               (1|survey_unit),
                             data = dissimilarity_covariates_dredge.dt)

global_mod_min_temp_lmer_surveyyear <- lmer(annual_dissimilarity_value ~ 
                             survey_unit:yearly_min_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit:summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                             season + yearly_mean_bypoint_avg.scaledacrossall + summed_tonnes_scaled_byreg +
                               (1|survey_unit) + (1|year),
                             data = dissimilarity_covariates_dredge.dt)

dissimilarity_covariates_dredge.dt[, survey_unit := factor(survey_unit, levels = sort(all_surveys))]

global_mod_min_temp_lmer_survey <-lmer(annual_dissimilarity_value ~ 
                             survey_unit:yearly_min_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit:summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                             season + yearly_mean_bypoint_avg.scaledacrossall + summed_tonnes_scaled_byreg +
                               (1|survey_unit),
                             data = dissimilarity_covariates_dredge.dt)

global_mod_seas_temp_lmer_surveyyear <- lmer(annual_dissimilarity_value ~ 
                             survey_unit:yearly_seas_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit:summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                             season + yearly_mean_bypoint_avg.scaledacrossall + summed_tonnes_scaled_byreg +
                               (1|survey_unit) + (1|year),
                             data = dissimilarity_covariates_dredge.dt)

global_mod_seas_temp_lmer_survey <- lmer(annual_dissimilarity_value ~ 
                             survey_unit:yearly_seas_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit:summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                             season + yearly_mean_bypoint_avg.scaledacrossall + summed_tonnes_scaled_byreg +
                               (1|survey_unit),
                             data = dissimilarity_covariates_dredge.dt)

global_mod_SD_temp_lmer_surveyyear <- lmer(annual_dissimilarity_value ~ 
                             survey_unit:yearly_mean_bypoint_SD.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit:summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                             season + yearly_mean_bypoint_avg.scaledacrossall + summed_tonnes_scaled_byreg +
                               (1|survey_unit) + (1|year),
                             data = dissimilarity_covariates_dredge.dt)

global_mod_SD_temp_lmer_survey <- lmer(annual_dissimilarity_value ~ 
                             survey_unit:yearly_mean_bypoint_SD.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit:summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                             season + yearly_mean_bypoint_avg.scaledacrossall + summed_tonnes_scaled_byreg +
                               (1|survey_unit),
                             data = dissimilarity_covariates_dredge.dt)

View(AICc(global_mod_mean_temp_lmer_survey, global_mod_max_temp_lmer_survey, global_mod_min_temp_lmer_survey, global_mod_seas_temp_lmer_survey, global_mod_SD_temp_lmer_survey,
          global_mod_mean_temp_lmer_surveyyear, global_mod_max_temp_lmer_surveyyear, global_mod_min_temp_lmer_surveyyear, global_mod_seas_temp_lmer_surveyyear, global_mod_SD_temp_lmer_surveyyear))

#build data table to report AICc
global_mod_temp_table <- data.table(
  `Random effect year` = c(
    c(rep(F,5),rep(T,5))
  ),
  `Temperature variable` = rep(c(                  
              "Average mean SBT",
              "Average maximum SBT",                        
              "Average minimum SBT",                  
              "Average SBT seasonality",                 
              "SBT SD"),2),
  deltaAICc = signif(min(AICc(global_mod_mean_temp_lmer_survey, global_mod_max_temp_lmer_survey, global_mod_min_temp_lmer_survey, global_mod_seas_temp_lmer_survey, global_mod_SD_temp_lmer_survey,
          global_mod_mean_temp_lmer_surveyyear, global_mod_max_temp_lmer_surveyyear, global_mod_min_temp_lmer_surveyyear, global_mod_seas_temp_lmer_surveyyear, global_mod_SD_temp_lmer_surveyyear)[,2])-AICc(global_mod_mean_temp_lmer_survey, global_mod_max_temp_lmer_survey, global_mod_min_temp_lmer_survey, global_mod_seas_temp_lmer_survey, global_mod_SD_temp_lmer_survey,
          global_mod_mean_temp_lmer_surveyyear, global_mod_max_temp_lmer_surveyyear, global_mod_min_temp_lmer_surveyyear, global_mod_seas_temp_lmer_surveyyear, global_mod_SD_temp_lmer_surveyyear)[,2],2))

#order by aicc
setorder(global_mod_temp_table,cols = -"deltaAICc")

global_mod_temp_table[,Rank := seq(1,10,by = 1)]

global_mod_sbt_table <- global_mod_temp_table[,.(Rank,`Random effect year`,`Temperature variable`, deltaAICc)]

fwrite(global_mod_sbt_table, here::here("output","global_mod_sbt_table.csv"))

Best performing global model includes minimum temperature (centered and scaled) (Still the case as of October 4, 2024)

Now, look at different combinations of all predictors (min temp) using dredge

Global model: global_mod_min_temp

colnames(dd.dt.2)

 [1] "(Intercept)"                                         "area_km.scaledacrossall"                            
 [3] "depth_annual_avg.scaledacrossall"                    "depth_annual_range.scaledacrossall"                 
 [5] "haul_id_count_annual.scaledacrossall"                "latitude_annual_avg.scaledacrossall"                
 [7] "latitude_annual_range.scaledacrossall"               "season"                                             
 [9] "spp_count_annual.scaledacrossall"                    "summed_tonnes_scaled_byreg"                         
[11] "yearly_mean_bypoint_avg.scaledacrossall"             "summed_tonnes_scaled_byreg:survey_unit"             
[13] "survey_unit:yearly_mean_bypoint_avg.scaledacrossall" "df"                                                 
[15] "logLik"                                              "AICc"                                               
[17] "delta"                                               "weight"                                             

Predict dissimilarity across years using averaged model (model_avg_delta4)

dissimilarity_covariates_dredge.dt_predictions <- copy(dissimilarity_covariates_dredge.dt)

  #allowing temp and fishing to vary (aka no changes)
  dissimilarity_covariates_dredge.dt_predictions[,pred_dissim := predict(model_avg_delta4, se.fit = T, full = F)[[1]]][,pred_se := predict(model_avg_delta4, se.fit = T, full = F)[[2]]] #full allows us to switch back to mixed effect models

Warning: argument 'full' ignoredWarning: argument 'full' ignored

  
  
  
  #constant temp in regions (aka take mean of temperature values within survey units so they are the same value within each year for a survey)
  dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg <- copy(dissimilarity_covariates_dredge.dt)
  dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg[,yearly_min_bypoint_avg.scaledacrossall:=mean(yearly_min_bypoint_avg.scaledacrossall),.(survey_unit)]

  #constant fishing in regions (aka take mean of fishing values within survey units so they are the same value within each year for a survey)
  dissimilarity_covariates_dredge.dt_predictions_consistentfishinginreg <- copy(dissimilarity_covariates_dredge.dt)
  dissimilarity_covariates_dredge.dt_predictions_consistentfishinginreg[,summed_tonnes_scaled_byreg:=mean(summed_tonnes_scaled_byreg),.(survey_unit)]

#and then with consistent temp and fishing in regions (aka take mean of both min temp and fishing)
dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg <- copy(dissimilarity_covariates_dredge.dt)
dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg[,summed_tonnes_scaled_byreg:=mean(summed_tonnes_scaled_byreg),.(survey_unit)][,yearly_min_bypoint_avg.scaledacrossall:=mean(yearly_min_bypoint_avg.scaledacrossall),.(survey_unit)]


#allowing temp and fishing to vary within regs (normal)
    dissimilarity_covariates_dredge.dt_predictions[,pred_dissim := 
                                                     predict(model_avg_delta4, se.fit = T, full = F, newdata =   dissimilarity_covariates_dredge.dt_predictions)[[1]]][,pred_se := 
                                                                                                                                                                         predict(model_avg_delta4, se.fit = T, full = F, newdata =   dissimilarity_covariates_dredge.dt_predictions)[[2]]]

Warning: argument 'full' ignoredWarning: argument 'full' ignored

  
#allowing only fishing to vary within regions (with mean temp)
  
   dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg[,pred_dissim := 
                                                                        predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg)[[1]]][,pred_se :=
                                                                                                                                                                                                              predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg)[[2]]]

Warning: argument 'full' ignoredWarning: argument 'full' ignored

  
  
#allowing only temperature to vary within regions (with mean fishing pressure)
     
   dissimilarity_covariates_dredge.dt_predictions_consistentfishinginreg[,pred_dissim :=
                                                                           predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistentfishinginreg)[[1]]][,pred_se := 
                                                                                                                                                                                                                                  predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistentfishinginreg)[[2]]]

Warning: argument 'full' ignoredWarning: argument 'full' ignored

   
#and then with consistent temp and fishing across regions (mean of both fishing and temp)
dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg <- copy(dissimilarity_covariates_dredge.dt)
dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg[,summed_tonnes_scaled_byreg:=mean(summed_tonnes_scaled_byreg),.(survey_unit)][,yearly_min_bypoint_avg.scaledacrossall:=mean(yearly_min_bypoint_avg.scaledacrossall),.(survey_unit)]

dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg[,pred_dissim := 
                                                                         predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg)[[1]]][,pred_se :=
                                                                                                                                                                                                                 predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg)[[2]]]

Warning: argument 'full' ignoredWarning: argument 'full' ignored

#and then with consistent temp and fishing across regions (mean of both fishing and temp)
dissimilarity_covariates_dredge.dt_predictions_consistenttempfishing <- copy(dissimilarity_covariates_dredge.dt)
dissimilarity_covariates_dredge.dt_predictions_consistenttempfishing[,summed_tonnes_scaled_byreg:=mean(summed_tonnes_scaled_byreg)][,yearly_min_bypoint_avg.scaledacrossall:=mean(yearly_min_bypoint_avg.scaledacrossall)]

dissimilarity_covariates_dredge.dt_predictions_consistenttempfishing[,pred_dissim := 
                                                                         predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistenttempfishing)[[1]]][,pred_se :=
                                                                                                                                                                                                                 predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistenttempfishing)[[2]]]

Warning: argument 'full' ignoredWarning: argument 'full' ignored

#FOR COLOR TO MATCH
#sort color link by survey name season
#alphabetical order
color_link_alpha <- setorder(color_link, survey_unit)

#exclude surveys we don't include
color_link_alpha <- color_link_alpha[survey_unit %in% unique(dissimilarity_covariates_dredge.dt_predictions$survey_unit),]

color_alpha_order <- color_link_alpha[,hex]
label_alpha_order <- color_link_alpha[,Survey_Name_Season]
  

#maintain temp and fishing
#plot
predicted_values_temp_fishing <- ggplot(dissimilarity_covariates_dredge.dt_predictions) +
  geom_point(aes(x = year, y = annual_dissimilarity_value, color = survey_unit)) +
  geom_line(aes(x = year, y = pred_dissim, color = survey_unit)) +
  geom_ribbon(aes(x = year, ymin = pred_dissim-pred_se, ymax = pred_dissim+pred_se, fill = survey_unit), alpha = 0.3) +
  scale_color_manual(values = color_alpha_order, labels = label_alpha_order, name = "Survey") +
  scale_fill_manual(values = color_alpha_order, labels = label_alpha_order, name = "Survey") +
  labs(x = "Year",y = "Average annual total\nBray Curtis dissimilarity") +
  ylim(0,1.5) +
  theme_classic() +
  ggtitle("Average model predictions")

#average temp and fishing for each region
predicted_values_temp_fishing_meantempfishinginsurvey <- ggplot(dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg) +
  geom_point(aes(x = year, y = annual_dissimilarity_value, color = survey_unit)) +
  geom_line(aes(x = year, y = pred_dissim, color = survey_unit)) +
  geom_ribbon(aes(x = year, ymin = pred_dissim-pred_se, ymax = pred_dissim+pred_se, fill = survey_unit), alpha = 0.3) +
  scale_color_manual(values = color_alpha_order, labels = label_alpha_order, name = "Survey") +
  scale_fill_manual(values = color_alpha_order, labels = label_alpha_order, name = "Survey") +
  labs(x = "Year",y = "Average annual total\nBray Curtis dissimilarity") +
  ylim(0,1.5) +
  theme_classic() +
  ggtitle("Average model predictions with mean\nsurvey temperature and fishing pressure")

#average temp and fishing across all regions
predicted_values_temp_fishing_meantempfishing <- ggplot(dissimilarity_covariates_dredge.dt_predictions_consistenttempfishing) +
  geom_point(aes(x = year, y = annual_dissimilarity_value, color = survey_unit)) +
  geom_line(aes(x = year, y = pred_dissim, color = survey_unit)) +
  geom_ribbon(aes(x = year, ymin = pred_dissim-pred_se, ymax = pred_dissim+pred_se, fill = survey_unit), alpha = 0.1) +
  scale_color_manual(values = color_alpha_order, labels = label_alpha_order, name = "Survey") +
  scale_fill_manual(values = color_alpha_order, labels = label_alpha_order, name = "Survey") +
  labs(x = "Year",y = "Average annual total\nBray Curtis dissimilarity") +
  ylim(0,1.5) +
  theme_classic() +
  ggtitle("Average model predictions with mean\novereall temperature and fishing pressure")

#merge plots
predicted_values_sbt_jaccard_fishing_merge <- cowplot::plot_grid(predicted_values_temp_fishing,
                                                 predicted_values_temp_fishing_meantempfishinginsurvey,
                                                 predicted_values_temp_fishing_meantempfishing,
                                                 ncol = 1)

Warning: Removed 27 rows containing missing values (`geom_line()`).Warning: no non-missing arguments to max; returning -InfWarning: no non-missing arguments to max; returning -Inf

ggsave(predicted_values_sbt_jaccard_fishing_merge, path = here::here("figures"), filename = "predicted_values_sbt_jaccard_fishing_merge.jpg", height = 30, width = 14)

Take dissimilarity values from random normal distribution for each year for each region, and then calculate slope (1000 times). Do this for: - Fishing and temperature vary interannually within surveys - Temperature is held constant (as mean over time series for a survey), but fishing varies (allows us to look at relative variance explained) - Fishing is held constant (as mean over time series for a survey), but temperature varies (allows us to look at relative variance explained) - Both fishing and temperature held constant (allows us to see role of other components of the model)

################################################################################
#full predictions
################################################################################
#table with predicted dissimilarity values and standard error of all predicted dissimilarity values (by year)
table <- dissimilarity_covariates_dredge.dt_predictions[,.(survey_unit, pred_dissim, pred_se, year)]
#0) make datatable to populate
  predicted_dissim_trends_rnormruns <- data.table()
#1) NEW PREDICTED VALUES FROM DISTRIBUTION
for (i in 1:1000){
  table[,rnorm_pred := rnorm(1, mean = pred_dissim, sd = pred_se),.(year, survey_unit)]
#2) CALCULATE LINEAR MODEL TO EXTRACT SURVEY SPECIFIC TREND VALUES
  jaccard_total_predicted_lm_singlerun <- lm(rnorm_pred ~ year*survey_unit,data = table)

  model_coefs_reduced_predictions_singlerun <- data.table(summary(jaccard_total_predicted_lm_singlerun)$coefficients)
  model_coefs_reduced_predictions_singlerun[,var := rownames(summary(jaccard_total_predicted_lm_singlerun)$coefficients)]
  
  #limit to interactions only (check this if there are any model changes!) row 2 and rows 34:64
  model_coefs_reduced_predictions_singlerun <- model_coefs_reduced_predictions_singlerun[c(2,34:64),]
  
  #adjust survey unit name by deleting beginning of string
  model_coefs_reduced_predictions_singlerun[,survey_unit := substr(var, 17, str_length(var))][var == "year",survey_unit := "AI"]
  
  #calculate interaction coefficients
  AI_estimate <- model_coefs_reduced_predictions_singlerun[1,Estimate]
  model_coefs_reduced_predictions_singlerun[1,estimate := AI_estimate]
  model_coefs_reduced_predictions_singlerun[2:32,estimate := (AI_estimate + Estimate)]
  
  predicted_dissim_trends_rnormruns <- rbind(predicted_dissim_trends_rnormruns, model_coefs_reduced_predictions_singlerun[,.(survey_unit, estimate)])
  
  print(i)
}
  
#reduce to mean and standard deviation
predicted_dissim_trends_rnormruns[,mean_dissim_coef:= mean(estimate),survey_unit][,sd_dissim := sd(estimate),.(survey_unit)]

predicted_dissim_trends_rnormruns.summary <- unique(predicted_dissim_trends_rnormruns[,.(survey_unit, mean_dissim_coef, sd_dissim)])

predicted_dissim_trends_rnormruns.summary[,pred_type := "full"]

################################################################################
#predictions with temp held constant and fishing still varying from year to year
################################################################################
#table with predicted dissimilarity values and standard error of all predicted dissimilarity values (by year)
table_constanttemp <- dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg[,.(survey_unit, pred_dissim, pred_se, year)]
#0) make datatable to populate
  predicted_dissim_trends_rnormruns_constanttemp <- data.table()
#1) NEW PREDICTED VALUES FROM DISTRIBUTION
for (i in 1:1000){
  table_constanttemp[,rnorm_pred := rnorm(1, mean = pred_dissim, sd = pred_se),.(year, survey_unit)]
#2) CALCULATE LINEAR MODEL FOR SLOPE VALUES
  jaccard_total_predicted_lm_singlerun_constanttemp <- lm(rnorm_pred ~ year*survey_unit,data = table_constanttemp)

  model_coefs_reduced_predictions_singlerun_constanttemp <- data.table(summary(jaccard_total_predicted_lm_singlerun_constanttemp)$coefficients)
  model_coefs_reduced_predictions_singlerun_constanttemp[,var := rownames(summary(jaccard_total_predicted_lm_singlerun_constanttemp)$coefficients)]
  
  #limit to interactions only (check this if there are any model changes!) row 2 and rows 34:64
  model_coefs_reduced_predictions_singlerun_constanttemp <- model_coefs_reduced_predictions_singlerun_constanttemp[c(2,34:64),]
  
  #adjust survey unit name by deleting beginning of string
  model_coefs_reduced_predictions_singlerun_constanttemp[,survey_unit := substr(var, 17, str_length(var))][var == "year",survey_unit := "AI"]
  
  #calculate interaction coefficients
  AI_estimate <- model_coefs_reduced_predictions_singlerun_constanttemp[1,Estimate]
  model_coefs_reduced_predictions_singlerun_constanttemp[1,estimate := AI_estimate]
  model_coefs_reduced_predictions_singlerun_constanttemp[2:32,estimate := (AI_estimate + Estimate)]
  
  predicted_dissim_trends_rnormruns_constanttemp <- rbind(predicted_dissim_trends_rnormruns_constanttemp, model_coefs_reduced_predictions_singlerun_constanttemp[,.(survey_unit, estimate)])
  
  print(i)
}
  
#reduce to mean and standard deviation
predicted_dissim_trends_rnormruns_constanttemp[,mean_dissim_coef:= mean(estimate),survey_unit][,sd_dissim := sd(estimate),.(survey_unit)]

predicted_dissim_trends_rnormrunsconstant_temp.summary <- unique(predicted_dissim_trends_rnormruns_constanttemp[,.(survey_unit, mean_dissim_coef, sd_dissim)])

predicted_dissim_trends_rnormrunsconstant_temp.summary[,pred_type := "temp_constant"]

predicted_dissim_trends_rnormruns.summary <- rbind(predicted_dissim_trends_rnormruns.summary, predicted_dissim_trends_rnormrunsconstant_temp.summary)

################################################################################
#predictions with fishing held constant (and temperature varying)
################################################################################
#table with predicted dissimilarity values and standard error of all predicted dissimilarity values (by year)
table_constantfishing <- dissimilarity_covariates_dredge.dt_predictions_consistentfishinginreg[,.(survey_unit, pred_dissim, pred_se, year)]
#0) make datatable to populate
  predicted_dissim_trends_rnormruns_constantfishing <- data.table()
#1) NEW PREDICTED VALUES FROM DISTRIBUTION
for (i in 1:1000){
  table_constantfishing[,rnorm_pred := rnorm(1, mean = pred_dissim, sd = pred_se),.(year, survey_unit)]
#2) CALCULATE LINEAR MODEL FOR SLOPE VALUES
  jaccard_total_predicted_lm_singlerun_constantfishing <- lm(rnorm_pred ~ year*survey_unit,data = table_constantfishing)

  model_coefs_reduced_predictions_singlerun_constantfishing <- data.table(summary(jaccard_total_predicted_lm_singlerun_constantfishing)$coefficients)
  model_coefs_reduced_predictions_singlerun_constantfishing[,var := rownames(summary(jaccard_total_predicted_lm_singlerun_constantfishing)$coefficients)]
  
  #limit to interactions only (check this if there are any model changes!) row 2 and rows 34:64
  model_coefs_reduced_predictions_singlerun_constantfishing <- model_coefs_reduced_predictions_singlerun_constantfishing[c(2,34:64),]
  
  #adjust survey unit name by deleting beginning of string
  model_coefs_reduced_predictions_singlerun_constantfishing[,survey_unit := substr(var, 17, str_length(var))][var == "year",survey_unit := "AI"]
  
  #calculate interaction coefficients
  AI_estimate <- model_coefs_reduced_predictions_singlerun_constantfishing[1,Estimate]
  model_coefs_reduced_predictions_singlerun_constantfishing[1,estimate := AI_estimate]
  model_coefs_reduced_predictions_singlerun_constantfishing[2:32,estimate := (AI_estimate + Estimate)]
  
  predicted_dissim_trends_rnormruns_constantfishing <- rbind(predicted_dissim_trends_rnormruns_constantfishing, model_coefs_reduced_predictions_singlerun_constantfishing[,.(survey_unit, estimate)])
  
  print(i)
}
  
#reduce to mean and standard deviation
predicted_dissim_trends_rnormruns_constantfishing[,mean_dissim_coef:= mean(estimate),survey_unit][,sd_dissim := sd(estimate),.(survey_unit)]

predicted_dissim_trends_rnormrunsconstant_fishing.summary <- unique(predicted_dissim_trends_rnormruns_constantfishing[,.(survey_unit, mean_dissim_coef, sd_dissim)])

predicted_dissim_trends_rnormrunsconstant_fishing.summary[,pred_type := "fishing_constant"]

predicted_dissim_trends_rnormruns.summary <- rbind(predicted_dissim_trends_rnormruns.summary, predicted_dissim_trends_rnormrunsconstant_fishing.summary)


################################################################################
#predictions with both fishing and temperature held constant (variability goes to other factors we don't account for)
################################################################################
#table with predicted dissimilarity values and standard error of all predicted dissimilarity values (by year)
table_constanttempfishing <- dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg[,.(survey_unit, pred_dissim, pred_se, year)]
#0) make datatable to populate
  predicted_dissim_trends_rnormruns_constanttempfishing <- data.table()
#1) NEW PREDICTED VALUES FROM DISTRIBUTION
for (i in 1:1000){
  table_constanttempfishing[,rnorm_pred := rnorm(1, mean = pred_dissim, sd = pred_se),.(year, survey_unit)]
#2) CALCULATE LINEAR MODEL FOR SLOPE VALUES
  jaccard_total_predicted_lm_singlerun_constanttempfishing <- lm(rnorm_pred ~ year*survey_unit,data = table_constanttempfishing)

  model_coefs_reduced_predictions_singlerun_constanttempfishing <- data.table(summary(jaccard_total_predicted_lm_singlerun_constanttempfishing)$coefficients)
  model_coefs_reduced_predictions_singlerun_constanttempfishing[,var := rownames(summary(jaccard_total_predicted_lm_singlerun_constanttempfishing)$coefficients)]
  
  #limit to interactions only (check this if there are any model changes!) row 2 and rows 34:64
  model_coefs_reduced_predictions_singlerun_constanttempfishing <- model_coefs_reduced_predictions_singlerun_constanttempfishing[c(2,34:64),]
  
  #adjust survey unit name by deleting beginning of string
  model_coefs_reduced_predictions_singlerun_constanttempfishing[,survey_unit := substr(var, 17, str_length(var))][var == "year",survey_unit := "AI"]
  
  #calculate interaction coefficients
  AI_estimate <- model_coefs_reduced_predictions_singlerun_constanttempfishing[1,Estimate]
  model_coefs_reduced_predictions_singlerun_constanttempfishing[1,estimate := AI_estimate]
  model_coefs_reduced_predictions_singlerun_constanttempfishing[2:32,estimate := (AI_estimate + Estimate)]
  
  predicted_dissim_trends_rnormruns_constanttempfishing <- rbind(predicted_dissim_trends_rnormruns_constanttempfishing, model_coefs_reduced_predictions_singlerun_constanttempfishing[,.(survey_unit, estimate)])
  
  print(i)
}
  
#reduce to mean and standard deviation
predicted_dissim_trends_rnormruns_constanttempfishing[,mean_dissim_coef:= mean(estimate),survey_unit][,sd_dissim := sd(estimate),.(survey_unit)]

predicted_dissim_trends_rnormrunsconstant_tempfishing.summary <- unique(predicted_dissim_trends_rnormruns_constanttempfishing[,.(survey_unit, mean_dissim_coef, sd_dissim)])

predicted_dissim_trends_rnormrunsconstant_tempfishing.summary[,pred_type := "fishing_and_temp_constant"]

predicted_dissim_trends_rnormruns.summary <- rbind(predicted_dissim_trends_rnormruns.summary, predicted_dissim_trends_rnormrunsconstant_tempfishing.summary)

Plotting observed vs predicted

jaccard_fishing_temp_model_observed_predicted_dt <- jaccard_total_coefs.r[predicted_dissim_trends_rnormruns.summary, on = "survey_unit"]

jaccard_fishing_temp_model_observed_predicted_dt[,pred_lower := mean_dissim_coef-sd_dissim][,pred_upper := mean_dissim_coef+sd_dissim]

#FULL MODEL, both temperature and fishing are allowed to vary
jaccard_fishing_temp_model_observed_predicted_lm <- lm(estimate ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "full"])
summary(jaccard_fishing_temp_model_observed_predicted_lm) #R^2 0.56

Call:
lm(formula = estimate ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == 
    "full"])

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0021560 -0.0005673 -0.0001628  0.0002739  0.0026908 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      4.504e-05  2.087e-04   0.216    0.831    
mean_dissim_coef 1.339e+00  2.140e-01   6.254 6.86e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.001179 on 30 degrees of freedom
Multiple R-squared:  0.566, Adjusted R-squared:  0.5515 
F-statistic: 39.12 on 1 and 30 DF,  p-value: 6.863e-07

(jaccard_fishing_temp_model_observed_predicted <- ggplot(jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "full"]) +
  geom_errorbar(aes(x = mean_dissim_coef, ymin = lwr, ymax = upr), color = "lightgrey", linewidth = 0.4) +
  geom_errorbarh(aes(y = estimate, xmin = mean_dissim_coef-sd_dissim, xmax = mean_dissim_coef+sd_dissim), color = "lightgrey", linewidth = 0.4) +
  geom_point(aes(y = estimate, x = mean_dissim_coef)) +
  geom_smooth(aes(y = estimate, x = mean_dissim_coef), color = "darkgrey",linetype = "dotted", method = "lm") +
  geom_abline(aes(slope = 1, intercept = 0)) +
  lims(x = c(min(jaccard_fishing_temp_model_observed_predicted_dt$pred_lower),max(jaccard_fishing_temp_model_observed_predicted_dt$pred_upper))) +
  labs(y = "Observed β-diversity trend",x = "Predicted β-diversity trend\n") +
  theme_classic()
)

#fishing constant (fishing constant; temperature varies only)
jaccard_fishing_constant_model_observed_predicted_lm <- lm(estimate ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "fishing_constant"])
summary(jaccard_fishing_constant_model_observed_predicted_lm) 

Call:
lm(formula = estimate ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == 
    "fishing_constant"])

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0038730 -0.0007399  0.0001136  0.0011372  0.0025014 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -0.0001537  0.0002541  -0.605 0.549774    
mean_dissim_coef  1.7942377  0.4072227   4.406 0.000124 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.001394 on 30 degrees of freedom
Multiple R-squared:  0.3929,    Adjusted R-squared:  0.3726 
F-statistic: 19.41 on 1 and 30 DF,  p-value: 0.0001238

#Temperature as a predictor, not fishing = R^2 = 0.40 (drop in 16% of variance explained when you lose fishing as predictor)

(jaccard_fishing_constant_model_observed_predicted <- ggplot(jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "fishing_constant"]) +
geom_errorbar(aes(x = mean_dissim_coef, ymin = lwr, ymax = upr), color = "lightgrey", linewidth = 0.4) +
  geom_errorbarh(aes(y = estimate, xmin = mean_dissim_coef-sd_dissim, xmax = mean_dissim_coef+sd_dissim), color = "lightgrey", linewidth = 0.4) +
  geom_point(aes(y = estimate, x = mean_dissim_coef)) +
    geom_smooth(aes(y = estimate, x = mean_dissim_coef), color = "darkgrey",linetype = "dotted", method = "lm") +
  geom_abline(aes(slope = 1, intercept = 0)) +
      lims(x = c(min(jaccard_fishing_temp_model_observed_predicted_dt$pred_lower),max(jaccard_fishing_temp_model_observed_predicted_dt$pred_upper))) +
  labs(y = "Observed β-diversity trend",x = "Predicted β-diversity trend\n(temperature varies fishing constant)") +
  theme_classic()
)

#temp constant (fishing only; temperature constant)
jaccard_temperature_constant_model_observed_predicted_lm <- lm(estimate ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "temp_constant"])
summary(jaccard_temperature_constant_model_observed_predicted_lm) #0.26 Drop in 30% of variance explained when you lose temperature as a predictor

Call:
lm(formula = estimate ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == 
    "temp_constant"])

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0025579 -0.0007197 -0.0002048  0.0003897  0.0042665 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)   
(Intercept)      0.0001684  0.0002729   0.617  0.54192   
mean_dissim_coef 0.8098204  0.2508463   3.228  0.00301 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.001542 on 30 degrees of freedom
Multiple R-squared:  0.2578,    Adjusted R-squared:  0.2331 
F-statistic: 10.42 on 1 and 30 DF,  p-value: 0.00301

(jaccard_temperature_constant_model_observed_predicted <- ggplot(jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "temp_constant"]) +
  geom_errorbar(aes(x = mean_dissim_coef, ymin = lwr, ymax = upr), color = "lightgrey", linewidth = 0.4) +
  geom_errorbarh(aes(y = estimate, xmin = mean_dissim_coef-sd_dissim, xmax = mean_dissim_coef+sd_dissim), color = "lightgrey", linewidth = 0.4) +
  geom_point(aes(y = estimate, x = mean_dissim_coef)) +
    geom_smooth(aes(y = estimate, x = mean_dissim_coef), color = "darkgrey",linetype = "dotted", method = "lm") +
  geom_abline(aes(slope = 1, intercept = 0)) +
      lims(x = c(min(jaccard_fishing_temp_model_observed_predicted_dt$pred_lower),max(jaccard_fishing_temp_model_observed_predicted_dt$pred_upper))) +
  labs(y = "Observed β-diversity trend",x = "Predicted β-diversity trend\n(fishing varies temperature constant)") +
  theme_classic()
)

#both temperature and fish held constant
jaccard_fishing_temp_model_observed_predicted_tempfishconstantinsurvey_lm <- lm(estimate ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "fishing_and_temp_constant"])
summary(jaccard_fishing_temp_model_observed_predicted_tempfishconstantinsurvey_lm) #%20 #drop of 36 from full

Call:
lm(formula = estimate ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == 
    "fishing_and_temp_constant"])

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0027797 -0.0008842 -0.0002533  0.0008447  0.0037051 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)  
(Intercept)      5.513e-05  2.846e-04   0.194   0.8477  
mean_dissim_coef 1.700e+00  6.284e-01   2.706   0.0111 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.001604 on 30 degrees of freedom
Multiple R-squared:  0.1962,    Adjusted R-squared:  0.1694 
F-statistic: 7.322 on 1 and 30 DF,  p-value: 0.01113

(jaccard_fishing_temp_model_observed_predicted_tempfishconstantinsurvey <- ggplot(jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "fishing_and_temp_constant"]) +
geom_errorbar(aes(x = mean_dissim_coef, ymin = lwr, ymax = upr), color = "lightgrey", linewidth = 0.4) +
  geom_errorbarh(aes(y = estimate, xmin = mean_dissim_coef-sd_dissim, xmax = mean_dissim_coef+sd_dissim), color = "lightgrey", linewidth = 0.4) +
  geom_point(aes(y = estimate, x = mean_dissim_coef)) +
    geom_smooth(aes(y = estimate, x = mean_dissim_coef), color = "darkgrey",linetype = "dotted", method = "lm")+
  geom_abline(aes(slope = 1, intercept = 0)) +
      lims(x = c(min(jaccard_fishing_temp_model_observed_predicted_dt$pred_lower),max(jaccard_fishing_temp_model_observed_predicted_dt$pred_upper))) +
  labs(y = "Observed β-diversity trend",x = "Predicted β-diversity trend\n(fishing and temperature constant)") +
  theme_classic()
)


#merge
jaccard_fishing_sbt_model_observed_predicted_merge <- plot_grid(jaccard_fishing_temp_model_observed_predicted + theme(plot.margin = unit(c(0.1,0.3,0.1,0.1),"cm")),
                                                                     jaccard_fishing_constant_model_observed_predicted + theme(plot.margin = unit(c(0.1,0.3,0.1,0.1),"cm")),
                                                                     jaccard_temperature_constant_model_observed_predicted + theme(plot.margin = unit(c(0.1,0.3,0.1,0.1),"cm")),
                                                                     jaccard_fishing_temp_model_observed_predicted_tempfishconstantinsurvey + theme(plot.margin = unit(c(0.1,0.3,0.1,0.1),"cm")), ncol = 2, labels = c("a.","b.","c.","d."))

`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'

ggsave(jaccard_fishing_sbt_model_observed_predicted_merge, path = here::here("figures"),filename = "jaccard_fishing_sbt_model_observed_predicted_merge.jpg", height =6, width = 8)

NA
NA

Let’s visualize model coefficients with temperature and fishing (similar to figure 2) model_avg_values

Plot all other coefficients in averaged model

model_avg_values.nonfishortemp <- model_avg_values[c(2,6,102:106),]

#mark significance
model_avg_values.nonfishortemp[,Significant := ifelse((coef-se > 0 & coef+se > 0) | (coef-se < 0 & coef+se < 0),T,F)]

#more helpful names for variables
model_avg_values.nonfishortemp[,Variable := c("Area","Species count","Depth","Number of tows","Latitude","Depth range","Latitude range")]

#make factor with order
model_avg_values.nonfishortemp[,Variable := factor(Variable, levels = c("Area","Species count","Number of tows","Depth","Depth range","Latitude","Latitude range"))]


#plot
all_avg_model_coef <- ggplot() + 
geom_point(data = model_avg_values.nonfishortemp, aes(x = Variable, y = coef, color = Significant)) +  
geom_errorbar(data = model_avg_values.nonfishortemp, aes(x = Variable, ymin = coef-se, ymax = coef+se, color = Significant), width = 0) +
  scale_color_manual(values = c("darkgrey","black")) +
geom_hline(yintercept = 0) +
scale_x_discrete(limits = rev) +
  labs(y = "Coefficient", x = "\n\n\n") +
coord_flip() +
theme_classic()

#plot
all_but_lat_avg_model_coef <- ggplot() + 
geom_point(data = model_avg_values.nonfishortemp[Variable == "Latitude"], aes(x = Variable, y = coef, color = Significant)) +  
geom_errorbar(data = model_avg_values.nonfishortemp[Variable == "Latitude"], aes(x = Variable, ymin = coef-se, ymax = coef+se, color = Significant), width = 0) +
  scale_color_manual(values = c("darkgrey","black")) +
geom_hline(yintercept = 0) +
scale_x_discrete(limits = rev) +
  labs(y = "Coefficient", x = "Model variable") +
coord_flip() +
theme_classic()

#plot
lat_avg_model_coef <- ggplot() + 
geom_point(data = model_avg_values.nonfishortemp[Variable != "Latitude"], aes(x = Variable, y = coef, color = Significant)) +  
geom_errorbar(data = model_avg_values.nonfishortemp[Variable != "Latitude"], aes(x = Variable, ymin = coef-se, ymax = coef+se, color = Significant), width = 0) +
  scale_color_manual(values = c("darkgrey","black")) +
geom_hline(yintercept = 0) +
scale_x_discrete(limits = rev) +
  labs(y = "Coefficient", x = "Model variable") +
coord_flip() +
theme_classic()

#merge into single plot

model_coef_summary_sbt_jaccard <- cowplot::plot_grid(sbt_fishing_avg_model_coef+theme(legend.position = "null", axis.title.x = element_blank()),all_avg_model_coef+theme(legend.position = "null"), ncol = 1, labels = c("  a.                                                      b.","        c."), label_y = 0.99, rel_heights = c(3,1))

ggsave(model_coef_summary_sbt_jaccard, path = here::here("figures"),filename = "model_coef_summary_sbt_jaccard.jpg", height = 6.5, width = 8, unit = "in")

#IF WE ACTUALLY USE THIS, NEED TO UPDATE ##Show intercepts by region and season # #{r} # ##survey intercepts #survey_intercepts <- model_avg_values[c(1,8:38),.(coef_name, coef, se)] #survey_intercepts[1,survey_unit := "AI"][1,coef_true := coef] #survey_intercepts[2:32,survey_unit := substr(coef_name, 12, str_length(coef_name))][2:32,coef_true := survey_intercepts[1,coef_true]+coef] #survey_intercepts <- color_link[survey_intercepts, on = "survey_unit"] #survey_intercepts[,Survey_Name_Season := reorder(Survey_Name_Season, coef_true)] # # ##season_intercepts #season_intercepts <- model_avg_values[c(1,3:5),.(coef_name, coef, se)] #season_intercepts[1,season := "Spring"][1,coef_true := coef] #season_intercepts[2:4,season := substr(coef_name, 7, str_length(coef_name))][2:4,coef_true := season_intercepts[1,coef_true]+coef] #season_intercepts[,season := reorder(season, coef_true)] # # ##survey intercepts #survey_model_coef <- ggplot() + #geom_point(data = survey_intercepts, aes(x = Survey_Name_Season, y = coef_true)) + #geom_errorbar(data = survey_intercepts, aes(x = Survey_Name_Season, ymin = coef_true-se, ymax = coef_true+se), width = 0) + #scale_x_discrete(limits = rev) + ## ylim(0.35,1) + # labs(y = "Intercept", x = "") + #coord_flip() + #theme_classic() # ##season intercepts #season_model_coef <- ggplot() + #geom_point(data = season_intercepts, aes(x = season, y = coef_true)) + #geom_errorbar(data = season_intercepts, aes(x = season, ymin = coef_true-se, ymax = coef_true+se), width = 0) + #scale_x_discrete(limits = rev) + # # ylim(0.35,1) + # labs(y = "Intercept", x = "") + #coord_flip() + #theme_classic() # ##merge into single plot # #model_intercept_jaccard <- cowplot::plot_grid(survey_model_coef+theme(axis.title.x = element_blank()), # season_model_coef, ncol = 1, labels = c("a.","b."), label_y = 0.99, rel_heights = c(3,1), align = "v") # #ggsave(model_intercept_jaccard, path = here::here("figures"),filename = "model_intercept_jaccard.jpg", height = 6.5, width = 6.5, unit = "in") # # # ##